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ABSTRACT 


This  study  compares  three  single  variable  search  methods  — 
Golden  Section,  cubic  interpolation  and  quadratic  interpolation. 
The  SUMT  nonlinear  program  was  used  for  the  comparison.  The 
OPT  subroutine  which  performs  the  single  variable  search  in 
SUMT  currently  uses  the  Golden  Section  method.  Two  different 
OPT  subroutines  were  written  which  implemented  cubic  inter¬ 
polation  and  quadratic  interpolation.  Seven  test  problems 
which  contained  9-100  variables  and  2-20  constraints  were  used. 
The  comparison  was  made  on  computation  time  per  single  variable 
search  for  the  three  methods  and  the  number  of  function  evalua¬ 
tions  per  single  variable  search  for  the  Golden  Section  and 
quadratic  interpolation  methods. 

A  single  variable  search  by  Lasdon,  Fox  and  Ratner  and 
one  by  Fletcher  and  McCann  were  also  discussed. 

The  results  showed  that  the  quadratic  interpolation  was 
slightly  faster  than  the  other  two  methods  and  required  fewer 
function  evaluations  per  single  variable  search  than  the 
Golden  Section  method.  Time  per  single  variable  search  was 
approximately  the  same  for  the  cubic  interpolation  and  Golden 
Section  methods.  Cubic  interpolation  required  fewer  points  to 
be  evaluated  than  the  other  two  methods  but  the  need  for 
gradient  evaluations  proved  to  be  costly  in  terms  of  computation 
time  per  single  variable  search. 
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I. 


INTRODUCTION 


The  concept  of  optimization  has  grown  to  play  a  major 
role  in  the  analysis  of  many  complex  decision  and  allocation 
problems.  The  quality  of  the  analysis  not  only  depends  upon 
the  skill  and  good  judgement  used  in  interpreting  and 
modeling  the  problem  but  also  upon  the  reliability  and 
efficiency  of  the  vehicle  used  for  finding  a  solution  to  the 
problem.  The  solving  of  nonlinear  programming  problems  has 
seen  substantial  development  during  the  past  twenty  years 
and,  no  doubt,  will  increase  in  acceptance  and  popularity  in 
the  future . 

Within  the  Department  of  Defense  several  problems  of 
significant  importance,  such  as  weapons  targeting  and  alloca¬ 
tion,  aircraft  and  power  plant  design  and  manpower  planning, 
to  name  a  few,  have  arisen  which  are  nonlinear  in  nature. 

Some  examples  of  industrial  nonlinear  programming  problems 
include  oil  refining,  curve  fitting,  inventory  and  logistics. 

Many  methods  exist  for  solving  nonlinear  programming 
problems.  Some  of  these  methods,  called  penalty  function 
methods,  are  based  on  transforming  a  constrained  minimization 
problem  into  a  sequence  of  unconstrained  minimization 
problems  whose  solutions  approach  the  solution  of  the 
original  constrained  problem. 

In  the  penalty  function  method,  as  well  as  in  most  other 
nonlinear  programming  algorithms,  a  one  dimensional 


8 


minimization  search  is  used  repeatedly  in  the  solution 
process.  In  this  paper  several  different  one  dimensional 
minimization  search  methods  are  compared  in  the  context  of 
penalty  function  algorithms . 

The  different  one  dimensional  minimization  search  methods 
(also  called  single  variable  searches)  are  compared  by  the 
amount  of  computation  time  and  by  the  number  of  fxinction, 
gradient  and  Hessian  evaluations  of  the  objective  function 
and  constraints  needed  to  find  a  solution  to  the  original 
constrained  optimization  problem. 

Section  II  of  this  paper  explains  the  formulation  of 
the  penalty  function  method  and  where  the  single  variable 
search  is  used  in  this  method.  The  nonlinear  computer 
program  SUMT  is  used  in  comparing ‘ the  different  single 
variable  search  methods  and  is  also  explained. 

Section  III  formulates  the  different  single  variable 
search  methods  —  Golden  Section,  ciibic  interpolation  and 
quadratic  interpolation  —  and  also  discusses  several  other 
methods  which  were  not  programmed  because  of  their 
complexity . 

Results  and  comparisons  of  the  different  methods  are 
stated  in  Section  IV  with  conclusions  and  suggestions  for 
further  study  in  Section  V.  The  Appendix  contains  the 
computer  codes  for  the  three  methods  that  were  programmed. 
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II.  THE  PENALTY  FUNCTION  METHOD 


A.  THE  OPTIMIZATION  PROBLEM 

The  single  variable  search  is  a  major  part  of  finding 
an  optimal  solution  to  a  constrained  optimization  problem 
by  the  penalty  function  method.  Dr.  W.C.  My lander,  an 
author  of  the  SUMT  computer  code,  estimated  that  forty  percent 
of  the  computation  time  required  to  find  an  optimization 
problem  solution  using  SUMT  involved  single  variable  searches. 

The  constrained  optimization  problem  to  be  solved  is  of 
the  form; 

minimize;  f(x) 

subject  to;  Uj (x)  ^0  j  =  l,...,m 

hj (x)  =  0  j  =  m+l,...,m+p  , 

where  x  is  an  n-dimensional  vector,  the  functions  f,  g  and  h 
are  continuous  and  may  be  either  linear  or  nonlinear,  and 
the  set  of  values  satisfying  the  inequality  constraints  has 
a  non-empty  interior,  i.e.  {x;  gj (x)  >  0,  j=l,...,m}  . 

B.  THE  SUMT  COMPUTER  PROGRAM 

The  computer  program  used  to  compare  the  different  single 
variable  search  methods  was  SUMT  -  Version  4  (Sequential 
Unconstrained  Minimization  Technique  for  Nonlinear  Programming) 
coded  in  FORTRAN  IV  by  W.C.  Mylander,  R.L.  Holmes  and  G.P. 
McCormick  for  the  Research  Analysis  Corporation  [Ref.  1] . 


10 


It  was  written  to  experiment  with  different  methods  of 
solving  nonlinear  programming  problems . 

A  Guide  to  SUMT  -  Version  4  [Ref.  2]  is  the  program  manual 
which  contains  detailed  information  about  the  development  of 
the  code,  the  subroutines,  options,  input-output,  limitations, 
user-supplied  subroutines,  data  deck  setup  and  a  detailed 
example  of  the  use  of  the  program  to  solve  a  nonlinear 
programming  problem. 

The  method  used  by  SUMT  to  solve  the  original  constrained 
optimization  problem  is  described  in  detail  in  Chapter  8  of 
the  book  on  nonlinear  programming  by  Fiacco  and  McCormick 
[Ref.  3] .  Basically,  SUMT  solves  the  constrained  problem 
by  finding  the  solutions  of  a  sequence  of  unconstrained 
problems  which  approach  the  solution  of  the  original  con¬ 
strained  problem. 

Previous  versions  of  SUMT  used  different  penalty  fTonctions 
for  approximating  the  solution  to  the  constrained  problem  but 
the  function  currently  used  to  sequentially  solve  the  problem 
is  of  the  form: 

m  m+p  2 

P(x,r)  =  f{x)  -  r  Z  In  g .  (x)  +  Z  [h.  {x)]^/r  , 

j=l  j=m+l 

where  the  parameter  r  determines  the  severity  of  the  penalty 
and  consequently  how  closely  the  unconstrained  problem 
solution  approaches  the  solution  to  the  constrained  problem. 

The  penalty  function  of  SUMT  uses  the  mixed  interior 
point-exterior  point  method  described  in  detail  in  Chapter  4 
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m 

of  Ref.  2.  The  term,  r  E  In  g • (x)  ,  involving  the 

,  j=l  ^ 

inequality  constraints  is  from  the  interior  point  method 

which  is  also  known  as  the  barrier  method.  As  the  x  vector 

approaches  the  boundary  of  the  infeasible  region  for  these 

constraints,  this  term  approaches  infinity.  The  term, 
m+P  2 

E  [h . (x) ]  /r  ,  involving  the  equality  constraints  is 
j=m+l  ^ 

from  the  exterior  point  method.  As  the  x  vector  moves 
farther  away  from  the  feasible  region  for  these  constraints, 
this  term  approaches  infinity.  The  behavior  of  these  two 
terms  is  shown  graphically  in  Figure  1.  How  much  of  a 
penalty  these  terms  add  to  the  penalty  function  depends  upon 
the  value  of  the  parameter  r. 

The  solution  procedure  requires  minimization  of  P(x,r) 
over  those  x  satisfying  the  conditions  gj (x)  >  0  ,  j=l,...,m 
for  r  =  r^,r2,...  where  r^^  <  r2  <...  <  rj^  <  ...  <0  .  Under 
mild  conditions  on  the  original  NLP  problem  the  minima  of 
the  sequential  unconstrained  problems  approach  the  solution 
of  the  original  constrained  problem  as  rj^  ^  0.  The  condi¬ 
tions  to  be  met  for  the  method  to  solve  the  original  problem 
are  described  in  detail  in  the  SUMT  program  manual  [Ref.  2]. 

When  the  NLP  problem  is  convex  the  SUMT  program  also 
generates  a  sequence  of  points  that  are  feasible  for  an 
optimization  problem  called  the  dual  of  the  original  prob¬ 
lem.  The.  maximum  function  value  of  the  dual  feasible  points 
and  minimum  function  value  of  the  unconstrained  problem 
bracket  the  optimal  solution  of  the  constrained  problem. 
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a.  INTEEIOR  POINT  METHOD  INVOLVING  E^EQUALITY  CCNSTRAniTS 


b.  EXTERIOR  POINT  METHOD  INVOLVING  EQUALITY  CCS'JSTRAINTS 
FIGURE  1.  BEHAVIOR  OF  THE  PENALTY  TERJ4S 
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As  the  unconstrained  problem  is  minimized  and  the  number  of 
dual  feasible  points  found  increases  the  difference  between 
the  two  solutions  decreases  and  the  optimal  solution  is  more 
accurately  approximated . 

SUMT  makes  use  of  the  theory  derived  in  Fiacco  and 
McCormick's  book  [Ref.  3,  Chapter  5]  that  uses  the  Lagrange 
extrapolation  technique  to  approximate  the  solution  of  the 
constrained  problem  when  more  than  one  minimum  of  the 
unconstrained  problem  has  been  found.  By  using  these 
extrapolation  approximations  the  solution  converges  faster 
and  achieves  a  closer  approximation  to  the  actual  solution 
than  the  minima  of  the  unconstrained  problems. 

The  logic  of  the  computer  program  taken  from  the  SUMT 
manual  [Ref.  2]  follows. 

STEP  1. 

Find  an  initial  starting  point  x^  within  the  interior 
feasible  region  such  that  gj (x)  >  0  ,  j=l,...,m  .  If  such 

an  initial  point  is  not  given  or  is  not  feasible/  the  program 
uses  the  optimization  method  to  find  one. 

j 

STEP  2. 

Determine  r^,  the  initial  value  of  r,  which  can  either 
be  an  input  parameter  or  a  rule  for  choosing  r^^  can  be 
specified. 

STEP  3. 

Determine  the  minimum  of  the  unconstrained  penalty 
function  for  the  current  value  of  r . 
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STEP  4. 


Estimate  a  better  solution  using  the  extrapolation 
technique,  if  possible. 

STEP  5. 

Computation  is  terminated  if  the  stopping  criteria  are 
satisfied  thus  yielding  an  approximate  solution  to  the 
o^iginsl  problem  with  accuracy  depending  upon  the  stringency 
of  the  stopping  criteria. 

STEP  6. 

If  the  stopping  criteria  are  not  satisfied  select 

^k+1  <  • 

STEP  7. 

Go  to  Step  3. 

A  simplified  flow  diagram  of  the  computer  logic  taken 
from  the  SUMT  manual  [Ref.  2,  p.  5]  is  shown  in  Figure  2. 

The  single  variable  search  problem  is  involved  in  Step  3 
of  the  procedure.  In  order  to  minimize  the  unconstrained 
penalty  function  the  program  chooses  a  search  direction  in 
n-space.  It  is  believed  that  by  searching  in  this  direction 
from  the  starting  point  of  the  subproblem  the  penalty  function 
will  initially  decrease.  A  one  dimensional  search  is  then 
made  along  that  direction  until  a  local  minimum  is  found. 

At  the  solution  to  this  one  dimensional  search  problem  a 
new  search  direction  is  computed  and  a  new  one  dimensional 
search  is  performed.  The  one  dimensional  searches  are 
repeated,  each  time  with  a  new  search  direction,  lantil  the 
unconstrained  problem  for  a  given  r  is  solved. 
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FIGURE  2. 


SUMT  Program  Flow  Diagram  [Ref.  2,  p.  5] 


SUMT  allows  the  user  to  choose  one  of  three  procedures 
to  determine  the  search  direction  —  by  Newton's  Method  using 
first  and  second  partial  derivatives  of  the  unconstrained 
penalty  function,  by  using  the  negative  of  the  gradient  of 
the  penalty  function  or  by  a  variable  metric  method  which  is 
taken  from  an  algorithm  of  the  Fiacco  and  McCormick  book 
[Ref.  3,  pp.  170-175]. 

Finding  a  solution  to  the  unconstrained  n  dimensional 
problem  thus  requires  several  single  variable  searches. 
Therefore,  the  more  efficient  the  single  variable  search  is, 
the  faster  the  computation  time  will  be.  Given  a  direction 
of  search,  the  SUMT  program  uses  the  subroutine  OPT  to 
perform  the  single  variable  search.  In  this  paper,  the 
presently  used  single  variable  search  method  in  the  OPT 
subroutine  is  compared  to  other  single  variable  search  methods . 
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III.  SINGLE  VARIABLE  SEARCH  METHODS 


Since  a  major  portion  of  finding  a  solution  to  the 
constrained  optimization  problem  involves  single  variable 
searches,  the  SUMT  program  was  used  to  compare  several 
single  variable  search  methods . 

The  current  OPT  subroutine  of  SUMT  uses  the  Golden 
Section  search  method  to  perform  the  single  variable  search. 
Two  different  single  variable  search  methods  were  programmed 
to  replace  the  current  OPT  subroutine.  The  first  method 

uses  cubic  interpolation  once  the  local  minimum  is 
bracketed  to  accelerate  convergence  to  that  local  minimum 
and  the  second  method  programmed  uses  quadratic  interpolation 
to  accelerate  convergence  once  three  feasible  points  that 
bracket  the  local  minimum  are  found.  Single  variable  search 
methods  by  Lasdon,  Fox  and  Ratner  [Ref.  4]  and  Fletcher  and 
McCann  [Ref.  5]  are  also  discussed. 

The  object  of  the  single  variable  search  is  to  find  a 
scalar  6  called  the  step  length  such  that  x*  =  +  0s 

where  x*  is  the  x  vector  that  corresponds  to  the  local 
minimum  of  the  penalty  function  along  the  search  direction  s 
from  an  initial  starting  point  x^.  For  ease  of  notation 
penalty  function  values  are  shown  as  a  function  of  step 
length  0  in  each  single  variable  search  since  the  initial 
starting  point  x^,  the  search  direction  s  and  the  parameter 
r  are  constant,  i.e.  P(0)  =  P(x^+0*s,r). 
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The  optimal  step  length  could  possibly  be  found  by 
moving  along  the  search  direction  in  a  random  manner. 
However,  this  would  not  be  a  very  logical  way  to  attack  the 
problem  and  would  lead  to  many  needless  evaluations  of  the 
penalty  function  which  would  increase  computation  time. 

The  OPT  subroutine  is  supplied  with  a  search  direction 
and  an  initial  point  from  which  to  search.  It  must  call 
other  subroutines  of  the  SUMT  program  for  functional, 
gradient  or  Hessian  evaluations  and  is  expected  to  return 
with  a  penalty  function  value  and  corresponding  x  vector 
which  is  the  local  minimum  for  that  single  variable  search. 
It  is  assumed  that  the  penalty  function  has  a  finite  local 
minimum  along  the  given  search  direction. 

The  following  subsections  formulate  and  explain  the 
different  single  variable  search  methods. 

A.  GOLDEN  SECTION  SEARCH  METHOD 

The  Golden  Section  search  method  uses  the  limiting 
properties  of  the  Fibonacci  series  [Ref.  6] .  Fiacco  and 
McCormick's  book  [Ref.  3]  define  the  steps  of  the  procedure 
used  in  the  Golden  Section  method. 

STEP  1. 

First  an  upper  bound  step  length  0^  is  obtained  with  the 
lower  bound  step  length  9^^  initially  equal  to  zero.  6^  is 
obtained  by  evaluating  the  penalty  function  at  successive 
step  lengths  which  are  in  the  limiting  Fibonacci  ratio, 
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(l  +  /5”)/2  ~  1.618  ,  that  is ,  6„  =  Z  (1.618)^  ,  where 

i=0 

k  is  the  smallest  integer  j  such  that 

3  0  o 

P  [  Z  (1.618)  ]  >  P  [  Z  (1.618) 

Jl=0  Jl=0 

with  6_  updated  as  Z  (1.618)  . 

^  )l=0 

STEP  2. 

The  interval  bounding  9*,  the  optimal  step  length/  is 
reduced  by  computing  two  other  step  lengths  within  the 
interval  (9^/9^).  Their  corresponding  values  are: 

9a  =  Bl  +  0.382  (9^  -  9^^) 

9j^  =  9j^  +  0 . 618  (  9y  -  9j^) 

STEP  3. 

The  penalty  function  values  of  the  two  interior  step 
lengths,  9^  and  9j^,  are  then  compared. 

STEP  4. 

If  P(9  )  <  P(9,),  then  the  optimal  step  length  which 
3.  iD 

corresponds  to  the  local  minimxam  of  the  single  variable 
search  should  be  in  the  interval  if  the  penalty 

function  is  unimodal.  Because  of  the  fact  that 
0.382/0.618  «  0.618,  reassign  9^  =  9j^,  9j^  =  9^  and  calculate 
a  new  9_  as  computed  in  Step  2.  Return  to  Step  3. 

3 
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STEP  5. 


If  ^(9^)  ^  P{6j^)f  then  the  optimal  step  size  should  be 

in  the  interval  (9^,0-.).  Reassign  0^  =  6  ,  9  =9^  and 

a.  u  j-i  a  a.  D 

calculate  a  new  0j^  as  computed  in  Step  2.  Return  to  Step  3. 
STEP  6. 

If  the  penalty  function  values  at  both  interior  step 

lengths  are  equal,  reassign  9,.  =  6  ,  0„  =  0,  and  calculate 

a  u  d 

a  new  0^  and  0j^  as  computed  before  by  returning  to  Step  2. 
STEP  7. 

When  the  stopping  criteria  is  satisfied  as  explained 
in  subsection  III.D,  0*  is  approximated  by  (0j^  +  0y)/2 
yielding  an  x  vector  x*  =  Xq+  9*S  which  corresponds  to  an 
approximated  local  minimum  P(0*)  of  the  single  variable 
search. 

A  flow  diagram  of  the  Golden  Section  method  presently 
used  in  the  OPT  subroutine  is  shown  in  Figure  3.  Certain 
parts  of  the  computer  code  that  were  the  same  in  each  of 
the  OPT  subroutines  are  discussed  in  subsection  III.D. 

The  Golden  Section  OPT  subroutine  initially  updates  0^ 

Lt 

as  it  finds  more  feasible  points  by  increasing  the  step 

length  until  it  finds  a  point  which  is  infeasible  or  where 

the  penalty  function  value  is  larger  than  the  penalty 

function  value  at  the  previous  feasible  point.  If  it  finds 

an  infeasible  step  length,  the  penalty  function  is  assigned 

3  6 

a  very  high  value  (10  ) ,  the  next  to  the  last  feasible 

step  length  is  assigned  0j^  and  the  last  feasible  step  length 
assigned  as  the  lower  interior  step  length  0  .  If  the  first 

a 
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FIGURE  3.  Golden  Section  OPT  Subroutine 
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step  length  from  the  initial  point  is  infeasible,  the  lower 

interior  step  length  9^  =  0.382  6^,  since  only  one  feasible 

step  length  (0  =  0)  has  been  found.  OPT  then  uses  0,,  and 

ij  U 

9^  to  compute  the  upper  interior  step  length 
9j2  =  9^  +  0 . 382  ( 9y  -  9^)  .  The  optimal  step  length  is 
bracketed  and  the  subroutine  continues  to  reduce  the 
interval  of  uncertainty  by  comparing  penalty  function  values 
of  the  interior  step  lengths  xintil  it  finds  two  values  or 
an  interval  that  satisfies  the  stopping  criteria. 

In  finding  the  minimum  of  each  single  variable  search, 
the  Golden  Section  method  only  requires  penalty  function 
evaluations  which  are  compared  to  show  where  the  penalty 
function  decreases  to  a  local  minimum  and  then  increases  to 
infinity  as  the  infeasible  region  boundary  is  approached. 

The  factor  by  which  the  interval  of  uncertainty  is  reduced 
remains  constant  in  this  method.  Therefore,  any  information 
about  the  behavior  of  the  penalty  function  other  than  the 
fact  that  the  optimal  step  length  is  somewhere  within  the 
interval  of  uncertainty  is  disregarded. 

B.  CUBIC  INTERPOLATION  METHOD 

The  c\abic  interpolation  method  was  programmed  for 
comparison  with  the  Golden  Section  and  quadratic  interpola¬ 
tion  methods.  Basically,  this  method  approximates  the 
penalty  function  by  a  cubic  fit  and  then  finds  the  minimum 
of  the  cubic  equation  to  approximate  the  optimal  step 
length  of  the  single  variable  search. 
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In  order  to  approximate  the  penalty  function  by  a  cubic 
fit  in  the  single  variable  search,  two  feasible  step  lengths 
which  bracket  the  local  minimum  along  with  their  directional 
derivative  values  are  needed.  This  method  actually  consists 
of  two  distinct  parts  —  the  bracketing  section  to  find  the 
two  feasible  points  and  the  cubic  interpolation  section  to 
find  a  minimiim  of  the  cubic  function.  The  single  variable 
search  method  is  supplied  with  an  initial  starting  point 
the  penalty  function  value  at  x^  [P(0)],  the  gradient  of 
the  penalty  function  at  x^  [VP(0)]  and  a  search  direction  s. 

The  steps  of  the  cubic  interpolation  method  follow. 

STEP  1. 

Since  the  directional  derivative  at  the  starting  point 
P' (0)  is  negative,  by  finding  a  step  length  that  has  a 
corresponding  positive  directional  derivative,  the  local 
minimum  along  the  given  direction  of  search  will  be  bracketed 
if  the  penalty  function  is  unimodal. 

The  step  length  is  increased  to  find  an  upper  step 
length  0^  until  an  n  is  found  such  that 

^  i 

P*  [(  Z  2^)  -  1]  >  0  . 

i=l 

If,  while  initially  increasing  the  step  length,  an  infeasible 
point  is  encountered,  the  increment  by  which  the  last  step 
length  was  increased  is  halved  and  subtracted  from  the 
current  step  length.  This  is  continued  until  a  step  length 
that  is  feasible  is  found.  If  its  directional  derivative  is 
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negative,  the  increment  is  increased  by  one  half  the 
increment  length  and  added  to  the  current  step  length. 

This  is  continued  until  a  step  length  with  a  positive 
directional  derivative  is  found. 

The  lower  step  length  6^  is  always  the  last  feasible 

Li 

step  length  encountered  with  a  negative  directional 
derivative . 

STEP  2. 

Once  two  step  lengths  have  been  found  that  bracket  the 
local  minimum  along  the  search  direction,  the  cubic  equation 
is  used  to  compute  a  step  length  which  corresponds  to  the 
minimum  of  the  cvibic  fit.  This  step  length  is  found  by  the 
equation: 


p'(e„)  +  Uj  •  °i 

e  -  -  (Sy  e^)  [  p' (e^)  -  +  20^ 


where 


P(9  -  POy) 

=  P'(e^)  H-  P'(e^j)  -  3  [  -  9^ . 


Uo  =  [u, 


-  p 


(0l)p 


(eg)] 


1/2 


STEP  3. 

P'(9)  is  evaluated  and  if  it  is  negative  then  reassign 
9,  =  8  and  the  left  side  is  discarded.  If  P'(6)  is  positive 

Li 

then  reassign  9^  =6  and  the  right  side  is  discarded. 
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STEP  4. 


If  the  stopping  criteria  is  satisfied,  the  two  penalty 
function  values,  P(9)  and  P(0j^)  or  P(e)  and  P  ( 0^)  ,  depending 
upon  P' (9) ,  are  compared  and  the  smallest  one  is  returned 
as  the  local  minimiam  P(0*)  of  the  single  variable  search 
with  its  corresponding  x  vector,  Xq  +  9*s. 

STEP  5. 

If  the  stopping  criteria  is  not  satisfied  return  to 
Step  2. 

A  flow  diagram  of  the  cubic  interpolation  OPT  subroutine 
is  shown  in  Figvire  4.  The  parts  of  this  subroutine  that  are 
the  same  in  the  other  OPT  subroutines  are  discussed  in 
subsection  III.D. 

The  ciabic  interpolation  method  uses  more  information  at 
each  feasible  step  length  to  find  the  local  minimvim  of  the 
single  variable  search  than  the  Golden  Section  or  quadratic 
interpolation  methods.  It  requires  a  gradient  evaluation 
for  the  penalty  function  at  each  feasible  step  length  in 
order  to  compute  the  directional  derivative.  A  gradient 
evaluation  may  require  many  function  evaluations  depending 
upon  the  number  of  variables  and  constraints  in  the  original 
constrained  problem.  However,  with  this  extra  information 
the  local  minimum  should  be  found  by  evaluating  fewer  step 
lengths  in  the  cubic  interpolation  section  of  this  method. 

As  the  value  of  the  parameter  r  decreases  with  each 
unconstrained  problem,  the  penalty  function  rises  more 
abruptly  as  the  infeasible  region  boundary  is  approached. 
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FIGURE  4.  Cubic  Interpolation  OPT  Subroutine 
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This  makes  it  more  difficult  to  find  a  step  length  with  a 
positive  directional  derivative  which  is  required  to  bracket 
the  minimum. 

C.  QUADRATIC  INTERPOLATION  METHOD 

A  quadratic  interpolation  single  variable  search  was  also 
programmed  for  comparison  with  the  two  previously  discussed 
methods . 

This  method  requires  three  feasible  step  lengths  and 

their  penalty  function  values  which  obey  the  relationship 

P(9^)  >  P(0^)  <  P(9,J  where  9,.  <  9,,  <  9„  .  With  this 
L  M  U  L  M  U 

information  a  quadratic  fit  is  made  with  its  minimum 
approximating  the  minimum  of  the  penalty  function  for  the 
given  search  direction.  The  penalty  function  is  evaluated 
at  the  interpolated  step  length  which  minimizes  the  quadratic 
function  and  compared  to  P  ( 9^^)  allowing  the  interval  of 
uncertainty  to  be  decreased.  The  process  is  repeated  until 
the  minimum  of  the  quadratic  fit  approximates  the  local 
minimum  of  the  penalty  function  or  the  interval  of  uncertainty 
becomes  small  enough  to  satisfy  the  stopping  criteria. 

This  method  also  consists  of  two  sections  —  the  bracketing 
section  and  the  quadratic  interpolation  section.  The  single 
variable  search  is  supplied  with  an  initial  x  vector  x^,  its 
penalty  function  value  P(0)  and  a  search  direction  s. 

The  steps  of  the  quadratic  interpolation  search  method 
follow. 
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STEP  1. 


The  step  length  is  increased  from  the  initial  starting 
point  to  find  an  upper  step  length  6^  until  an  n  is  found 
such  that 

n  .  n-1 

P[  (  Z  2^)  -  1]  >  P[  (  Z  2^)  -  1]  . 

i=l  i=l 

If,  while  initially  increasing  the  step  length,  an  infeasible 
step  length  is  encountered,  the  increment  by  which  the  last 
step  length  was  increased  is  halved  and  subtracted  from  the 
current  step  length.  This  is  continued  until  a  feasible  step 
length  is  found.  If  its  fiinctional  value'  is  less  than  P(9j^), 
the  increment  is  increased  by  one  half  the  increment  leng-th. 
This  is  continued  iintil  a  feasible  step  length  is  found  that 
has  a  larger  penalty  function  value  than  the  previous  step 
length.  Before  this  is  repeated,  reassign  0^^  =  0j^  and 

ensuring  that  * 

If  only  one  feasible  step  leng-th  is  found  in  obtaining 
the  upper  step  leng-th,  where  in  this  case  P(9u)  ^ 

(which  indicates  the  minimum  is  bracketed) ,  the  last 
increment  is  halved  and  subtracted  from  the  upper  step 
leng-th  which  yields  the  interior  step  length  0^. 

STEP  2. 

Now  that  the  local  minimxim  is  bracketed,  -the  quadratic 
function  is  used  to  find  a  step  length  where  the  derivative 
of  -the  quadratic  fit  vanishes. 
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i/j  —  1/2,3  where  1  =  L,  2  —  M,  and  3  =  U. 

STEP  3. 

If  0  >  0^^,  then  the  left  side  is  discarded  and  reassign 
9^  “  ®M  ~  If  0  <  then  the  right  side  is 

discarded  and  reassign  0,.  =  0..  and  0.  =  0. 

U  M  M 

STEP  4, 

P{0j^)  is  evaluated  and  if  the  stopping  criteria. are 
satisfied,  the  smallest  of  P{0j^),  P{0jyj)  and  P(0u)  is  returned 
as  P(0*)  with  the  corresponding  x  vector,  x^  +  0*s. 

STEP  5. 

If  the  stopping  criteria  are  not  met,  return  to  Step  2. 

A  flow  diagram  of  the  quadratic  interpolation  method  OPT 
subroutine  is  shown  in  Figure  5 . 

The  quadratic  interpolation  method  only  requires  function 
evaluations  to  perform  the  single  variable  search.  It  uses 
the  knowledge  of  three  feasible  step  lengths  and  their 
penalty  fvinction  values  to  approximate  where  the  local 
minimum  of  the  single  variable  search  is  located.  However, 
it  must  find  an  upper  feasible  step  length  that  is  within  the 
interval  where  the  penalty  fvmction  is  increasing  to  infinity 
and  also  this  feasible  step  length  must  have  a  penalty 
function  value  that  is  greater  than  that  of  the  interior 
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FIGUEE  5 .  Quadratic  Interpolation  OPT  Subroutine 
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feasible  step  length.  This  causes  the  interval  containing 
possible  feasible  upper  step  lengths  to  be  reduced  even  more. 
As  the  parameter  r  decreases  for  each  unconstrained  sub¬ 
problem,  the  penalty  function  increases  more  abruptly  from 
the  minimum  as  the  infeasible  region  boundary  is  approached. 
This  also  increases  the  difficulty  in  bracketing  the  minimum. 

Once  the  three  feasible  step  lengths  are  found,  the 
quadratic  interpolation  can  converge  to  the  local  minimum 
thereby  solving  the  single  variable  search. 

D.  SECTIONS  COMMON  TO  THE  DIFFERENT  METHODS 

Some  parts  of  the  three  different  OPT  subroutines  were 
made  similar  so  that  they  could  be  compared  under  the  same 
conditions. 

Essentially  the  same  stopping  criteria  was  used  in  each 
subroutine.  When  the  ratios  of  the  two  x  vectors  that 
bracketed  the  minimum  or  ratios  of  their  penalty  function 

-7 

values  were  within  10  of  one,  the  subroutine  returned  to 
the  program  with  a  local  minimum  and  a  corresponding  x  vector 
as  the  solution  to  the  single  variable  search.  These  criteria 
were  checked  in  the  Golden  Section  method  whenever  two 
interior  step  lengths  had  been  computed  and  in  the  inter¬ 
polation  methods  both  after  the  minimum  had  been  bracketed 
and  also  after  each  interpolated  step  length  had  been 
computed . 

In  the  three  different  methods,  the  subroutines  also 
contained  counters  which  would  cause  the  search  uo  stop 
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after  ten  interpolations  in  the  cubic  interpolation  method, 
twenty  interpolations  in  the  quadratic  interpolation  method 
or  twenty-five  feasible  points  in  the  Golden  Section  method. 
In  the  test  problems  that  were  run  these  counters  were  never 
reached,  thereby,  never  satisfying  this  stopping  criterion. 

In  the  cubic  interpolation  method  another  stopping 
criterion  was  used  along  with  the  ratio  tests.  When  the 
cosine  of  the  angle  between  the  gradient  of  the  penalty 

-7 

function  and  the  search  direction  vector  was  less  than  10 
the  subroutine  returned  with  a  local  minimum  and  x  vector. 

It  was  stated  earlier  that  it  was  hoped  that  the  penalty 
function  would  initially  decrease  in  the  given  search 
direction  from  the  initial  starting  point.  In  the  three 
different  methods,  if  the  ratio  of  any  element  of  the 
current  x  vector  and  the  corresponding  element  in  the  initial 

-7 

X  vector  came  within  10  of  one,  the  single  variable  search 
was  restarted  in  the  opposite  search  direction.  In  the 
interpolation  methods  if  after  one  hundred  tries  the  minimum 
was  not  bracketed  or  the  step  length  had  been  decreased  more 
than  twenty  consecutive  times,  the  search  direction  was  also 
reversed  and  the  single  variable  search  was  restarted. 

The  three  different  OPT  subroutines  called  the  subroutine 
EVALU  for  penalty  function  evaluations .  EVALU  would  return 
back  a  variable  after  each  call  up  which  would  show  if  the 
step  length  was  feasible.  The  cubic  interpolation  method 
called  the  subroutine  GRAD  for  the  gradient  of  the  penalty 
function  when  it  was  needed  to  determine  the  directional 
derivative. 
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E.  OTHER  SINGLE  VARIABLE  SEARCH  METHODS 

Two  other  single  variable  search  methods  were  also 
looked  at  but  were  not  programmed  because  of  their  complexity 
and  their  need  for  changing  other  portions  of  the  SUMT 
program  than  just  the  OPT  subroutine. 

1.  Lasdon^  Fox,  and  Ratner  Method 

The  single  variable  search  developed  by  Lasdon,  Fox 
and  Ratner  [Ref.  4]  uses  a  penalty  function  that  doesn't 
contain  equality  constraints  and  the  penalty  term  for  the 
inequalities  is  also  different  than  the  SUMT  penalty  function. 
The  unconstrained  minimization  problem  is: 

...  1 
minimize  P(x,r)  =  F(x)  +  r  Z  - — ?— y 

i=l 

where  F  is  the  objective  function  and  the  G^'s  are  the 
inequality  constraints.  This  model  could  be  modified  to 
handle  equality  constraints  and  a  penalty  function  like  that 
used  in  SUMT. 

The  single  variable  search  method  is  in  three  distinct 
stages  —  linear  approximation,  quadratic  approximation  and 
cubic  interpolation.  Linear  approximations  are  made  for  the 
objective  fianction  and  each  constraint  using  F(0),  F'(0), 

Gj(0)  and  Gj'(O).  A  step  length  (0^^)  which  minimizes  the 
penalty  function  consisting  of  the  linear  approximations  is 
then  calculated.  The  information  used  to  make  the  linear 
approximations  is  then  used  along  with  F(9j^),  F'(e^),  G^O^^) 
and  Gj'Oj^)  to  make  quadratic  approximations  of  the 
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objective  fxanction  and  constraints.  The  penalty  function 
consisting  of  quadratic  approximations  is  then  minimized 
by  computing  Q2^  If  the  directional  derivative  of  the 
original  penalty  function  is  positive  at  ©2/  cubic  inter¬ 
polation  is  applied  yielding  6*,  the  solution  to  the  single 
variable  search.  If  the  stopping  criterion  is  satisfied 
during  any  stage,  the  single  variable  search  is  terminated 
yielding  a  local  minimum  and  a  solution  to  the  search. 

The  following  steps  outline  the  Lasdon,  Fox  and 
Ratner  procedure. 

STAGE  1. 

The  objective  function  and  inequality  constraints 
of  the  original  constrained  problem  are  approximated  by 
using  the  given  values  at  the  initial  starting  point. 

f^(9)  =  F(0)  +  F’  (0) 

g^j(9)  =  Gj(0)  +  Gj'(O)  j  =  l,...,m. 

The  approximating  function  for  P(9)  is 

m 

=  f^(e)  +  r 

The  smallest  positive  zero  of  the  linear  approximations  of 
the  inequality  constraints  is  found  by  taking  the  smallest 
-Gj (0) /Gj ' (0)  over  all  j.  This  value  is  designated  the 
upper  step  length  9y. 

The  step  length  corresponding  to  the  minimxam  of  the 
approximated  penalty  function  is  found  by  doing  four  or  five 
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Newton's  method  iterations  using  ^  6^ 


as  the  starting  point. 


0 


1 


Pi"  <1  V 


The  result  of  this  stage  is  a  step  length  that 
approximates  the  optimal  step  length  of  the  single  variable 
search.  If  the  stopping  criterion,  explained  following 
Stage  3,  is  not  satisfied  proceed  to  Stage  2. 

STAGE  2. 

The  same  procedure  is  performed  in  this  stage  as  in 
Stage  1  with  the  exception  that  quadratic  approximations  to 
the  objective  function  and  constraints  are  made  instead  of 
linear  approximations.  The  second  stage  is  entered  with 
values  of  F(6j^)  and  Gj  (9^)  along  with  the  information  used 
to  make  the  linear  approximations.  The  quadratic  approxima¬ 
tions  to  the  objective  fiinction  and  constraints  are 


where 


f2(e)  =  ae^  +  be  +  c 


F(9^)  -  be^^  -  c 


b  =  F'  (0)  ,  c  =  F(0) 


and 


where 


92.(6)  =  d.e^  ^ 

G  (6^)  -  e.6i  - 

u  .  =  ^ - := - 

:  6,2 

e.  =  G.'(O)  ,  f.  =  G. (0). 
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If  Gj(9^)  is  infeasible  for  some  j,  then  the  smallest 
positive  root  over  all  j  of 


-e .  ±  /e -  4d.f . 

-J - 2_ 2 1 

2dj 


0 


2U 


is  used  as  the  starting  step  length  for  using  Newton's 
method  to  minimize  the  approximated  penalty  function 


m 

p,(9)  =  f, (0)  +  r  Z 
^  i=l 


The  step  length  which  minimizes  P2(0)  should  be 
found  after  four  or  five  iterations  of  Newton's  method. 


9 


2 


P2"<V 


If  is  infeasible,  use  j  02^  in  Newton's  method. 
If  the  stopping  criterion  is  not  satisfied  after  finding 
02  continue  to  Stage  3. 

STAGE  3. 

Direct  cubic  interpolation  of  the  penalty  function 
is  used  in  this  stage.  If  P'(02)  is  positive,  the  minimum 
is  bracketed  and  ciibic  interpolation  can  be  performed  to 
find  a  step  length  which  corresponds  to  the  minimum  of  the 
penalty  fxinction.  If  the  minimum  has  not  been  bracketed, 
that  is,  if  P'(92)  <  0,  a  bracketing  procedure  like  the  one 
used  in  the  cubic  interpolation  method  of  sxibsection  II I. B 
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must  be  performed  to  find  a  step  length  that  has  a  positive 
directional  derivative.  Once  this  step  length  is  found, 
cubic  interpolation  is  used  to  find  a  step  length  which 
minimizes  the  cubic  function.  If  the  new  step  length  does 
not  satisfy  the  stopping  criteria,  additional  cubic  inter¬ 
polations  are  made  using  the  two  points  which  most  closely 
bracket  the  minimum  of  the  penalty  function. 

When  two  penalty  function,  values  are  nearly  equal 
or  when  the  interval  between  step  lengths  that  bracket  the 
minimum  becomes  very  small,  an  optimal  step  length  is  deter¬ 
mined  by 

p'  (e^)  -  p'  (e. ) 

a  D 

where  6  is  the  lower  step  length  and  0,  is  the  upper  step 
length  (0^  <  0*  <  0j^)  . 

The  stopping  criterion  used  in  this  method  is  the 

I  s'^VP  I 

— I - 1 -  ,  where  (p  is  the  angle 

l|s|ll|VP|| 

between  the  search  direction  vector  and  gradient  of  the 
penalty  function.  When  the  cosine  of  the  angle  is  less  than 

—  O 

10  ,  the  stopping  criterion  is  satisfied. 

The  authors  found  that  in  using  this  method  on  test 
problems,  the  stopping  criterion  was  often  satisfied  at  the 
end  of  Stage  2.  In  their  comparisons  of  this  method  with  a 
cubic  interpolation  method,  they  found  that  it  performed 
the  single  variable  search  much  more  efficiently  [Ref.  4, 


cosine  test  |  cos  (j)  [ 


This  method  requires  gradients  of  the  objective 
function  and  constraints.  In  the  SUMT  program  these  values 
are  not  stored  so  implementing  this  method  would  require 
more  storage  or  more  gradient  evaluations  along  with 
changes  in  other  subroutines  and/or  addition  of  new 
subroutines . 

2 .  Fletcher-McCann  Method 

The  Fletcher-McCann  method  [Ref.  5,  pp.  210-212] 
uses  an  approximation  of  the  original  unconstrained  penalty 
function  to  find  the  local  minimum  of  the  single  variable 
search.  The  penalty  frinction  approximation  used  is  of  the 
form 


T(0) 


a  +  b0  +  c9 


2 


9 


where  the  parameters  a,  b,  c,  d  and  0^  are  determined  by 
using  objective  function  and  constraint  information  at  two 
feasible  step  lengths.  The  authors  felt  that  an  approximation 
of  this  type  which  goes  to  infinity  at  the  barrier  0  =  0^ 
would  fit  the  penalty  function  better  than  a  polynomial 
approximation  which  goes  to  infinity  only  as  0 
an  estimate  of  the  intersection  of  the  search  direction  s 
with  the  boundary  of  the  infeasible  region.  The  approximated 
penalty  function's  minimum  is  found  by  applying  Newton's 
method  to  find  a  step  length  which  corresponds  to  the  local 
minimum  along  the  direction  of  search. 

Since  the  explanation  of  this  method  by  the  authors 
was  rather  sketchy,  certain  parts  could  be  interpreted 
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^iffsirently  and.  only  rssolved  by  testing  the  method  on 
dif fstent  problems .  An  explanation  of  the  notation  for 
this  method  is  first  given. 


1  m+P  .  , 

F{e)  =  f(e)  +  i  E  h/(e)  ::  a  +  be  + 

j=m+l  ^ 


m 

S 

j=l 


GO)  =  -r  Z  in  ,.0)  = 


F' (6)  = 


P  p  ra+p 

Vf  (6)  +  ^  E  h.  (e)  s^Vh.  (9) 
^  i=m+l  ^ J 

Moll 


=:  b  +  C0 


G' (6)  =  - 


m  s  Vg . ( 6) 
E  - 3- 


j=i  93(e) 


(Sy  -  6) 


The  single  variable  search  is  supplied  with  an  initial  point, 
f\anctional  and  gradient  information  at  that  point  and  a 
search  direction.  The  steps  of  the  method  follow. 

STEP  1. 

A  feasible  step  length  in  the  direction  of  search  is 
needed  along  with  the  initial  point  in  order  to  estimate  the 
parameters  a,  b,  c,  d  and  6y.  A  unit  step  length  is  initially 
taken  in  trying  to  find  a  feasible  step  length.  If  it  is  not 
feasible,  the  step  length  is  progressively  halved  until  a 
feasible  step  length  is  found.  The  lower  and  upper  feasible 
step  lengths  are  written  as  6^  and  0^^,  respectively.  At  the 
start  of  the  search  9q  =  0.  Function  and  gradient  evaluations 
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are  made  at  the  second  feasible  step  length  in  order  to 
compute  F^,  F^'  and  G^'  (shortened  notation  for  F(0j^), 

G(0^) ,  etc. )  . 

By  using  the  information  at  the  two  feasible  step 
lengths,  simultaneous  equations  are  solved  to  determine 
the  parameters .  The  parameter  values  are 


U 


=o'  -  Si 


The  parameter  0y  is  taken  as  the  smallest  value  of  two 
computed  values  which  is  greater  than  0^^. 

STEP  2. 

Once  the  parameters  are  determined,  the  equation  for 
the  minimum  of  T  which  is 

T*  =  0  =  b  +  2c0  +  - j 

{0U  -  0) 

is  solved  by  Newton's  method  using  the  midpoint  between  the 
two  feasible  step  lengths  as  the  starting  point. 
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e 


9,  -  0. 


] 


where  T"  =  2c  +  ^  ^  . 

^  (0u“®) 

STEP  3. 

The  values  of  VP{0)  are  computed  and  if  the  stopping 
criterion  is  satisfied,  the  single  variable  search  is 
terminated  with  0  as  the  optimal  step  length. 

STEP  4. 

If  the  stopping  criterion  is  not  satisfied,  a  new 
step  length  is  computed  using  Newton's  method  in  Step  2. 

The  authors  of  this  method  stated  that  different 
default  actions  would  be  necessary  in  cases  where  the  inter¬ 
polation  failed  or  was  \insatisf actory .  For  example,  if  for 
the  feasible  step  lengths,  0q  and  0^^,  a  value  of  0^^  could 
not  be  determined,  it  would  be  necessary  to  use  another 
feasible  step  length  to  determine  where  the  intersection  of 
the  search  direction  and  the  boundary  of  the  infeasible 
region  was.  In  later  stages  of  the  minimization  the  authors 
said  that  this  method  sometimes  failed  and  a  quadratic 
interpolation  would  have  to  be  used. 

Implementing  this  method  in  the  SUMT  program  would 
also  require  more  storage  or  more  gradient  evaluations 
along  with  modifications  to  other  subroutines  and/or  additional 
subroutines . 


42 


No  comparisons  are  made  of  the  Lasdon,  Fox  and 
Ratner  or  Fletcher-McCann  methods  to  the  other  three 
methods  since  they  are  only  discussed  and  not  programmed. 
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IV.  TEST  PROBLEMS 


A.  TEST  PROBLEM  ORIGIN 

The  problems  used  to  compare  the  three  different  single 
variable  search  methods  programmed  as  OPT  subroutines  were 
taken  from  problems  used  in  a  thesis  by  Lt.  J.  Waterman,  USN, 
which  compared  three  different  nonlinear  programming  codes 
[Ref.  7].  The  SUMT  program  and  problem  data  decks  were  the 
same  as  those  used  by  Waterman  except  for  the  Golden  Section 
OPT  subroutine  being  replaced  by  the  cubic  and  quadratic 
interpolation  methods . 

The  structure  and  degree  of  difficulty  among  the  seven 
test  problems  are  quite  varied  and  represent  a  sample  of 
some  real  world  problems.  The  number  of  variables  and 
constraints  ranged  from  9  to  100  and  2  to  20  respectively. 

The  problems  contain  combinations  of  linear,  nonlinear, 
equality  and  inequality  constraints . 

The  following  subsection  contains  the  descriptions  of 
the  test  problems  as  presented  in  Waterman's  thesis. 

Constants  in  each  problem  are  represented  by  a,  b,  c,  d,  e,  L, 
m,  u  and  s  unless  otherwise  specified. 

B.  TEST  PROBLEM  DESCRIPTION 

Problem  1  was  an  example  of  determining  the  chemical 
composition  of  a  complex  mixture  under  conditions  of 
chemical  equilibrium.  It  contained  45  independent  variables 
and  16  linear  equality  constraints. 
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minimize  f(x) 


+  In 


X 


k=l  j=l 


iJL 


^  ^-ik 
j=l 


-)  ] 


siabject  to 


h^(x) 


7 

E  (  Z  E.  .j^x  ,  ) 
k=l  j=l 


=  0, 


1  ■“ 


^jk  -  ^  ^  l,...,hj^  k  -  1, - ,1 . 


Problem  2  was  formulated  by  the  Shell  Development 
Company.  It  consisted  of  15  variables  and  5  nonlinear  equality 
constraints . 


10  5  5  53 

maximize  f(x)  =  E  b.x.  -  Z  E  yz  -  2  E  d.z 

i=l  j=l  i=l  j=l  ^ 


where  y  =  and  z  =  ==(io+j) 


10 


subject 


to  g.(x)  =2  E  y  +  3  d.z  +  e .  -  E  a. .x.  ^  0 
3  i=l  i=i  13  1 


x^^O  i—l|..ofl5 


Problem  3  was  to  maximize  the  area  of  a  hexagon  in  which 
the  maximum  diameter  was  unity.  The  problem  had  9  independent 
variables,  13  nonlinear  inequality  constraints  and  a  lower 
boxind  of  0  for  Xg . 
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Maximize:  f{x)  -  •  5 (x^^x^  -  X2X2  +  x^x^  -  x^x^  +  x^Xg  -  XgX^) 

subject  to: 


2  2 

1  -  Xg  -  x^  ^0 


1  -  X4  >0 

1  -  x_^  -  •X.?'  >  0 
3  6  — 

1  -  (Xj^  -  Xg)  ^  -  {X2.-Xg)^  >  0 

2  2 

1  -  (Xj^-X^)  -  (X2-Xg)  ^  0 

1  -  {Xj^  -  x^)  ^  -  {X2  -  Xg)  ^  ^  0 

1  -  {Xg  -  Xg)  ^  -  {X4  -  Xg)  ^  ^  0 

1  -  (Xg  -  X^)  ^  -  (X4  -  Xg)  ^  ^  0 

1  -  X^^  -  (Xg  -  Xg)  ^  ^  0 


Xj^X^  -  XgXg  0 


X3X9  >  0 


-XgXg  >  0 


x-Xq  -  x^x-,  >  0 

5  8  6  7  — 


Xg  >_  0  . 


Problem  4  was  probably  the  most  difficult  of  the  test 
problems.  It  included  a  linear  objective  function,  24 
variables,  12  nonlinear  equality,  2  linear  equality  and 
6  non-linear  inequality  constraints.  The  variables  had 
to  also  be  zero  or  positive. 
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minimize; 


24 

f (x)  =  E  a.  X. 

i=l  ^  ^ 


subject  to; 


h^(x)  = 

X 

(i+12) 

CiXi 

24  X. 

b  Z  — ^ 

o 

cr 

X  . 

_J. 

b. 

3 

24 
=  Z 
i=l 

X.  -  1  =  0 

1 

h^4{x) 

12 
=  Z 
i=l 

X.  24 

-^  +  f  1 

i=13 

^ 

1 

=  0 

where 

f  =  142.224 

U 

(X)  = 

-(Xi  +  ^(^1+12) 

) 

4-  ^  >0 

.  1 

^(i+14) 

24 

i  ^  i  ^ 

=  0,  i  =  1, ... ,12 


=  1,2,3 


j=l 


X  . 

3 


24 

Z  X. 

j=l  ^ 


+  e^  ^  0,  i  —  4,5,6 


x^  ^  0 ,  i  =  1, . . . , 24 
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Problem  5  was  a  weapon  assignment  problem  with  100 
variables,  a  nonlinear  objective  function,  12  linear 
constraints  and  zero  lower  bounds  for  the  variables. 


minimize ; 


f  (x) 


20 

E 

j=l 


u . 


5 

n 

i=l 


1) 


subject  to; 


5 

Z  X. .  -  b.  >  0  j  =  1,6,10,14,15,16,20 
i=l  3  ~ 


20 

-  E  ^0  i  =  l,...,5 

j=l 


i  “  j  —  1 


20  . 


Problem  6  was  adapted  from  an  inventory  model  where 
the  x^,  i  =  1,...,50,  represent  the  reorder  quantity  for 
50  inventory  items  and  x^,  i  =  51,..., 100,  represent  the 
reorder  points  for  the  same  50  items.  It  contained  100 
variables,  1  linear  and  1  nonlinear  inequality  constraint, 
and  50  lower  bounds  on  the  variables. 
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minimize; 


50  B.  (x.  +  50) 
f(x)  =  E  ^ 

i=l 


X. 

1 


where 


B.(x^  +  50)  =i  (s.2+d.2)  $(^)  - 


s.d.  d. 

1  1  A  /  lx 

~T-  ' 

1 


^(i+50)  ” 


i-  e""" 


''2tt 


and  $(r)  =  /  (t)(x)  dx  . 

t 


subject  to: 


50  X. 

200,000  =i<^+-,i^50) 


-  m^)  ^  0 


50  L. 

300  -  E  —  >  0 

.  ,  X.  — 

1=1  1 


Xj^  ^  0  ,  i  =  1,  .  .  . ,  50  . 


Problem  7  was  an  entropy  model.  There  were  46  population 
centers  connected  by  a  transportation  network.  Using  a 
congestion  cost  function  the  model  yields  an  equilibrium 
solution  that  identifies  nodal  populations  as  entropic 
functions  of  the  total  cost  of  the  journey  to  work.  The 
problem  contained  46  variables,  all  of  which  have  a  zero 
lower  bound,  1  nonlinear  inequality  and  1  linear  equality 
constraint. 
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f  (x) 


46 

Z 

i=l 


minimize: 


Xi 

500  5^  ^ 


46 

subject  to:  500  -  Z  x-  =  0 

i=l  ^ 


10000  -  Z  c.y.  +  ad.y.  >  0 
i=l  ^  ^  ^  ^  - 

where  y.  =  x.  +  Z  x.  ,  A(i)  consists  of  all  the 

jeA(i)  ^ 

arcs  that  converge  directly  and  indirectly  upon  node  i, 
and 

x^^O  i=l,...,46  . 

The  preceding  problem  descriptions  are  only  meant  to 
give  a  feeling  of  the  kind  of  test  problems  used  and  their 
structure.  Detailed  information  of  how  each  problem  was 
set  up,  the  constant  values,  initial  starting  points  and 
where  the  problems  specifically  originated  is  contained 
in  Ref.  5. 
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V.  TEST  RESULTS  AND  COMPARISONS 


A.  PROGRAMMING  AND  TESTING  PROCEDURE 

The  cvibic  interpolation  method  was  first  programmed 
as  a  program  that  did  a  single  variable  search  when  given 
a  penalty  function,  its  gradient,  an  initial  starting  point 
and  a  search  direction.  After  it  had  been  debugged  it  was 
converted  into  an  OPT  subroutine.  After  the  subroutine 
was  debugged  using  a  couple  of  the  less  difficult  test 
problems,  OPT  was  converted  into  the  quadratic  interpolation 
method  by  modifying  the  bracketing  procedure  and  changing 
the  interpolation  from  a  cubic  function  to  a  quadratic 
function. 

When  the  quadratic  interpolation  method  had  been  debugged 
each  problem  was  run  through  the  SUMT  program  three  times , 
using  the  three  different  search  methods.  As  discussed 
earlier,  the  stopping  criteria  was  made  essentially  the  same 
in  each  of  the  OPT  subroutines  so  that  the  only  thing  that 
varied  in  solving  each  test  problem  was  the  single  variable 
search  method. 

In  the  process  of  running  the  other  test  problems,  minor 
debugging  changes  had  to  be  made  in  the  interpolation  methods 
After  all  the  changes  had  been  incorporated,  a  final  run 
of  each  problem  with  each  method  was  made  which  produced  the 
results  used  in  the  comparisons.  Each  method  found  the  same 
solution  to  each  problem  to  within  six  significant  figures. 
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B.  COMPARISON  CRITERIA 


t 


Perhaps  the  best  way  to  compare  the  different  search 
methods  is  to  compare  the  computation  time  required  to  find 
the  solution  to  the  problem.  However,  because  of  computer 
interactions,  computation  time  may  vary  as  much  as  twenty 

> 

five  per  cent  when  the  same  problem  is  run  at  two  different 
times .  Computation  times  were  computed  for  each  problem  to 
see  if  a  trend  could  be  seen  between  the  methods. 

Counters  were  inserted  in  the  SUMT  program  so  that  the 
number  of  single  variable  searches  performed  and  the  number 
of  functional  and  gradient  evaluations  needed  to  find  the 
solution  could  be  counted. 

Results  of  the  test  problem  runs  are  shown  in  Table  I. 

In  running  problem  7  with  the  interpolation  methods  it  was 
necessary  to  change  the  factor  by  which  the  increment  was 
multiplied  for  step  length  increase  or  reduction  to  1.618 
vice  2  and  0.618  vice  0.5  respectively.  This  change  was 
needed  because  a  feasible  starting  point  was  found  in  the 
bracketing  procedures  which  the  SUMT  program  could  not  solve. 
By  changing  the  factors  to  the  same  as  those  in  the  Golden 
Section  method,  the  same  feasible  starting  point  was  used 
to  solve  the  problem. 

C.  COMPARISONS 

Theoretically,  each  single  variable  search  should  reach 
the  same  solution  for  all  three  methods.  However,  because 
of  the  tolerances  allowed  in  the  stopping  critera,  the 
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TABLE  I 


TEST  PROBLEM  RESULTS 


Problem 

Golden  Section 

Cubic 

Quadrat! 

TIME 

105.8 

105.7 

94.8 

T 

SVS 

86 

89 

75 

X 

F 

15487 

6171 

10149 

G 

1462 

7497 

1275 

TIME 

11.4 

9.6 

8.3 

o 

SVS 

63 

53 

53 

F 

4354 

1574. 

2606 

G 

378 

1704 

318 

TIME 

4.8 

5.6 

4.0 

SVS 

39 

38 

40 

3 

F 

8606 

3302 

5862 

G 

613 

3358 

628 

TIME 

395.3 

404.1 

396.1 

A 

SVS 

163 

151 

174 

4 

F 

44255 

13922 

30048 

G 

92.27 

97319 

98694 

TIME 

103.7 

97.9 

97.0 

SVS 

20 

19 

19 

5 

F 

2683 

1250 

1991 

G 

296 

1415 

271 

TI11E 

165.4 

145.9 

155.5 

SVS 

36 

34 

35 

6 

F 

2320 

836 

978 

G 

104 

817 

107 

TIME 

78.8 

80.3 

67.8 

SVS 

20 

20 

17 

7 

F 

865 

307 

397 

G 

1993 

2257 

1708 

TIME  =  Computation  time 
SVS  =  #  single  variable  searches 
F  =  #  function  evaluations 
G  =  #  gradient  evaluations 
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solutions  were  slightly  different  at  the  conclusion  of  each 
single  variable  search.  This  meant  that  at  the  start  of 
the  next  single  variable  search  the  search  direction  would 
be  slightly  different  also.  As  several  single  variable 
searches  were  performed  in  each  problem,  the  differences 
accumulated  and  this  explains  why  for  some  of  the  problems 
it  takes  a  different  number  of  searches  for  each  method  to 
reach  the  same  solution.  This  also  explains  the  difference 
in  the  number  of  function  and  gradient  evaluations  required 
for  each  problem  for  the  Golden  Section  and  quadratic 
interpolation  methods. 

Table  II  shows  the  normalized  results  for  the  time  per 
single  variable  search,  function  evaluations  per  single 
variable  search  and  gradient  evaluations  per  single  variable. 

As  was  expected,  the  nximber  of  gradient  evaluations  per 
single  variable  search  is  essentially  the  same  for  Golden 
Section  and  quadratic  interpolation -methods .  The  only  thing 
that  can  be  compared  between  the  three  methods  is  the 
computation  time  per  single  variable  search.  Fxanction 
evaluations  per  single  variable  search  can  only  be  used  as 
a  measure  of  effectiveness  for  Golden  Section  and  quadratic 
interpolation  since  the  cubic  interpolation  also  requires 
gradient  evaluations  and  fewer  fxanction  evaluations. 

In  looking  at  the  times  per  single  variable  search  the 
quadratic  interpolation  was  faster  than  the  Golden  Section 
and  cubic  interpolation  methods  5  out  of  7  and  4  out  of  7 
test  problems  respectively.  Cubic  interpolation  was  faster 
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Cm  O 


TABLE  II.  NORMALIZED  TEST  PROBLEM  RESULTS 


Problem 

Golden  Section 

Cubic 

Quadratic 

TIME 

1.23 

1.19 

1.26 

1 

F 

180.1 

69.3 

135.3 

G 

17 

84.2 

17 

TIME 

.181 

.181 

.157 

2 

F 

69.1 

29.7 

49.2 

G 

6 

32.1 

6 

TIME 

.123 

.147 

.1 

3 

F 

220.7 

86.9 

146.6 

G 

15.7 

88.4 

15.7 

TIME 

2.42 

2.67 

2.28 

4 

F 

271.5 

92.2 

172.7 

G 

565.2 

644.5 

567.2 

TIME 

5.18 

5.15 

5.10 

5 

F 

64.4 

24.6 

27.9 

G 

2.8 

24.0 

3.0 

TIME 

4.59 

4.29 

4.44 

6 

F 

64.4 

24.6 

27.9 

G 

2.8 

24.0 

3.0 

TIME 

3.94 

4.02 

3.98 

7 

F 

43.3 

15.4 

23.4 

G 

99.7 

112.9 

100.4 

TIME  =  time/single  variable  search 

=  #  function  evaluations/single  variable  search 
=  #  gradient  evaluations/single  variable  search 
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than  Golden  Section  in  3  of  the  7  problems  with  the  same 
time  for  problem  2 . 

In  comparing  the  number  of  function  evaluations  per 
single  variable  search,  the  quadratic  interpolation  method 
required  fewer  than  the  Golden  Section  on  all  seven  problems. 

In  most  of  the  test  problems,  it  required  between  thirty  to 
forty  per  cent  fewer  function  evaluations.  The  cvibic  inter¬ 
polation  method  required  fewer  function  evaluations  than  the 
other  two  methods  because  it  used  more  information  about  the 
penalty  function  at  each  feasible  step  length  to  find  the 
optimal  step  length.  However,  it  required  more  gradient 
evaluations  which  increased  the  computation  time  making  it 
about  the  same  as  the  Golden  Section  and  slightly  slower  than 
the  quadratic  interpolation. 

With  a  larger  sample  of  test  problems  more  statistically 
sound  comparisons  could  be  made  of  the  three  different  methods. 
Also  samples  of  test  problems  that  are  similar  in  structure 
could  be  tested  to  show  if  one  method  worked  better  than  the 
others  on  those  type  of  problems. 

It  should  be  emphasized  that  the  results  in  this  thesis 
are  obtained  from  single  variable  searches  on  a  very  special 
type  of  function  —  the  xmconstrained  penalty  function  in  the 
mixed  interior— exterior  penalty  function  algorithm.  No 
attempt  should  be  made  to  generalize  these  results  to  other 
nonlinear  programming  methods  which  also  use  single  variable 
searches,  but  on  a  significantly  different  class  of  functions. 
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VI.  CONCLUSIONS  AND  SUGGESTIONS  FOR  FURTHER  STUDY 


A.  CONCLUSIONS 

In  looking  at  the  computation  time  required  for  each 
single  variable  search,  the  quadratic  interpolation  method 
was  faster  than  the  other  two  methods  on  the  test  problems 
that  were  run.  The  quadratic  interpolation  method  uses 
information  about  the  penalty  function  at  each  feasible  step 
length  to  reduce  the  interval  of  uncertainty  whereas  the 
Golden  Section  reduces  the  interval  of  uncertainty  by  a 
constant  factor.  Because  of  this  difference,  quadratic 
interpolation  proved  to  be  slightly  more  efficient  than 
Golden  Section.  The  only  drawback  in  the  quadratic  inter¬ 
polation  method  is  the  fact  that  bracketing  the  minimum  can 
be  quite  hard  to  do.  The  cubic  interpolation  method  reduced 
the  interval  of  uncertainty  in  a  more  efficient  manner  than 
quadratic  interpolation  or  Golden  Section.  However,  having 
to  make  gradient  evaluations  to  determine  the  directional 
derivative  values  proved  to  be  very  costly.  The  time  per 
single  variable  search  was  about  the  same  as  the  Golden 
Section  search  method. 

If  the  user  of  the  single  variable  search  method  has 
prior  knowledge  about  the  complexity  of  the  function  or 
gradient  evaluations  he  may  prefer  cubic  interpolation  over 
Golden  Section  or  quadratic  interpolation,  or  vice  versa, 
since  the  cubic  interpolation  requires  more  gradient  evalua¬ 
tions  and  fewer  function  evaluations. 
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B.  SUGGESTIONS  FOR  FURTHER  STUDY 

A  combination  of  different  methods  made  into  a  single 
variable  search  could  be  done  with  the  three  different 
methods  that  were  programmed.  This  combined  methods  approach 
could  use  Golden  Section  until  the  minimum  was  bracketed, 
then  use  quadratic  interpolation  until  the  minimum 
was  bracketed  much  closer  and  then  finally  use  the  cubic 
interpolation  to  home  in  on  the  local  minimum  of  the  search. 
The  Lasdon,  Fox  and  Ratner  method  uses  three  different  stages 
that  each  find  a  feasible  step  length,  with  the  final  stage 
finding  an  optimal  step  length  for  the  search.  Another 
approximation  of  the  penalty  function  like  the  Fletcher- 
McCann  method  could  also  be  devised. 

Implementation  of  the  Lasdon,  Fox  and  Ratner  or  Fletcher- 
McCann  method  into  the  SUMT  program  with  runs  of  the  test 
problems  would  enable  a  comparison  to  be  made  not  only  to 
the  three  methods  tested  in  this  paper  but  also  between  the 
Lasdon,  Fox  and  Ratner  method  and  Fletcher-McCann  method. 

Another  possibility  for  study  would  be  to  fine  tune  the 
three  different  methods  by  adjusting  the  tolerances  and 
termination  conditions  to  see  that  if  by  allowing  a  less 
accurate  determination  of  the  minimum  of  the  single  variable 
search,  the  solution  to  the  unconstrained  problem  could  be 
found  any  faster. 

If  a  larger  sample  of  test  problems  including  more  and 
varied  types  of  nonlinear  programming  problems  could  be 
tested  using  the  three  different  methods ,  a  better  comparison 
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could  be  made.  The  results  may  either  show  one  method  to 
be  superior  over  the  other  methods  in  all  of  the  problems 
or  show  each  method  suited  best  to  a  specific  type  of 
problem  enabling  the  user  to  choose  the  method  which  best 
applies  to  the  situation  at  hand. 
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APPENDIX 


DELX 

DELXO 

DOTT 

I 

ISW 

J 

KSW 

M 

MN 

N 

NSATIS 

NTCTR 

N401 

N404 

N405 

PREV3 


GLOSSARY  FOR  GOLDEN  SECTION  OPT  SUBROUTINE  [Ref.  2] 


The  vector  indicating  the  direction  of  move 
in  one  dimensional  optimization.  S 

The  gradient  vector  of  the  penalty  function 
VP. 


The  inner  product  of  the  move  vector  and  the 
negative  grac^ent  vector  of  the  penalty 
function.  -sTVp 


Used  as  an  index  in  DO  loops. 

A  switch  showing  whether  the  motion  on  a  given 
vector  failed  to  decrease  the  penalty  function 
and  the  negative  was  tried. 

Used  as  an  index  in  DO  loops . 

Switch  showing  whether  less  than  25  feasible 
points  were  found  on  the  direction  vector. 

The  nximber  of  inequality  constraints 

Number  of  moves  in  search  for  a  solution  of 
a  subproblem 

The  number  of  variables 

Indicates  whether  constraints  are  satisfied. 

The  number  of  the  point  on  which  the  program 
is  working. 

The  number  of  points  generated  in  attempting 
to  find  a  point  along  the  search  direction. 

The  number  of  infeasible  points  found  on  the 
direction  vector. 

Switch  showing  that  a  feasible  point  on  the 
direction  vector  could  not  be  found  and  the 
negative  was  tried. 

Penalty  function  value  at  lower  step  length 
P(X3) . 
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PXl 


Penalty  function  value  for  left  interior  steplength 
P(X2). 

PO  Current  value  of  the  penalty  function  P(X) 

Pi  Penalty  function  value  for  upper  steplength 

P(X1) 

P31  Penalty  function  value  at  initial  starting  point 

RJ  Vector  of  current  values  of  the  constraints 

RJl  Vector  of  previous  value  of  the  constraints 

X  Vector  of  current  value  of  the  variables  Xq+0s 

XX  Temporary  storage  used  for  switching  vector  values 

XI  X  vector  of  upper  step  length  Xq+9^s 

X2  X  vector  of  left  interior  step  length 

X3  X  vector  of  lower  step  length 
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SLERCUTINE  OPT 

^  I^FLICIT  REAL*8( A-h,0-Z) 

C  NARCH  1971 

C 

r  ALONG  THE  SEARCH  VECTOR  USING  THE 

C  GCLOEN  SECTION  SEARCH  McTHQC 

CCHPGN/SHARE/  X(IOO),  CEL(100)t  A  ( 100, IOC  )  » N , M ,  HN,NF1 

1  »  ^  r  1 

CCNJJCN  /VALUE/  F ,  G  ,  PO  ,R  S IGMA  ,  RJ(200),  RFC 
CCFPGN/CRST/  DcLX(lOO),  DELXO(IOO),  RHOIN, RATIO,  EPSI, 
ITHETAO, 

2  FSIGl,  Gl,  XKIOO),  X2(100),  >3(100),  XF2(100), 
3XP1  (100  )  ,PP1, 

^  RJ1(200),  DOTT,  PGRAC(lOO),  CIAG(IOO), 

5  FPEV3,ACELX,  NTCTR,  NUNINI,  NPHASE,  NSA7IS 
AES(CUMNY)=DABS(CUMNY) 

KSV»  =  1 
N^C5=1 
F 31=P0 
ISW=1 
CCT7=0o 
CC  10  J=1,N 

CCT7  =  C.TT  +  DELX(  J)’i‘0£LX0(J  ) 

GC  70  40 
CC  30  1=1, N 
CELX( I ) =-OELX( I ) 

CONTINUE 
N AC 4=0 
RN=NN+1 

C  MN  IS  NOW  NUMB.  GF  POINTS  AFTER  MIN  ACHIEVED 
NTCTR=NTCTR+1 
CC  50  1=1, N 
50  X2(I)=X(I) 

FX1=P0 
N4C1=0 
N AC1=N4C1+1 
CC  70  1=1, N 
X(I  )  =  X2 (I )+D£LX( I ) 

CALL  EVALU 

1  HEANS  SATIS. GF  CONSTRAINT  NT.PREV.  2MEANS  NCCHANGE  3 
MEANS  VIOLATION 

IF  POINT  IS  NCT  FEASIBLE  GIVE  IT  AN  ARBRITRARILV  HIGH  VALU 
GC  TO  (540,90,80),  NSATIS 
F>2=10.E25 
FC=10.E35 
GC  TO  100 
CONTINUE 
F>2=P0 

IF  (FX1-PX2)  100,100,150 
IF  (N401-2)  130,110,110 
CC  120  1  =  1, N 


10 

20 

30 

40 


60 

70 

C 

C 

C 

60 


90 


100 

110 

120 


SO  FAR  COMPUTED 


XKI  )=X(I  ) 

F  1=PX2 
GC  TO  430 
C  ONLY  CNE  PCINT 
130  CC  140  1=1, N 
140  X3(I)=X2(I) 

FREV3=P>1 
GC  TC  160 

150  CC  160  1  =  1, N 
X2( I)=X2( I ) 

160  CELX(i  )  =1.61803399=»CELX(I  ) 
FFEV3  =  P  XI 
F)1=PX2 
GC  TO  60 

C  GCLDEN  SECTION  SEARCH  METHOD. 

C  E  VECTOR  GCES  TO  XI  (I  ) 

170  PC=1.E36 

^404=N404+l 
180  CC  190  1  =  1, N 


190  X1{I)=X(I) 

P1  =  F0 

CC  2C0  I=1»N 

X(I)=,38196601#(X1(I)-X3(  1) )+X3(I) 

200  X2(I»=X(I) 

CiSLL  EV/iLU 

GC  TG  (540,270,210),  NS/MIS 

210  IF  (N404.LT. 30)  GC  TQ  170 

211  CONTINUE 

C  ThERE  IS  NC  REFERENCE  TC  211,  ^H5  AECVE  STATEMENT  IS  A 
C  OLMMT  STATEMENT 

C— IT  IS  POSSIBLE  NO  FEASIBLE  POINT  EXIST,  IF  NCT  TRY  MOVING 
C  ON  CELXQ 

C—  IF  IT  IS  NOT  POSSIBLE  TG  MOVE  ON  CELXC  THEN  Vc  MUST  EE 
C  AT  A  SOLUTION  OF  NLP  PROBLEM. 

IF  (N404.GT.100)  GO  TO  240 
220  CC  230  1  =  1, N 

IF  (AeS(ABS(X3( I  )/Xl( I )  )-l,  ),GT.l.E-7)  GC  TC  170 
220  CONTINUE 
240  C-C  TC  (250,260),  N405 
250  N4C5=2 

C — 7PY  TC  MCVe  CN  GRADIENT. 

NTCTR=NTCTR-1 
MN=MN-1 
GC  TO  20 

260  WRITE  (6,580) 

CALL  TIMEC 
CALL  OUTPUT  (1) 

CALL  REJECT 
STCP  22042 
C 

270  CONTINUE 
N4C4=0 
FX1=P0 

CC  280  1  =  1, N 

280  X (I )=0. 28 196601* ( XI (I )-X2(I))+X2(I) 

CALL  EVALU 

GC  TO  (540,290,220),  NSATIS 


290  FX2=P0 
N4C1=1 

200  N4C1=N401+1 

IF  (N401-25)  340,210,210 
210  KSW=2 

IF  (N4C1-40)  320,460,460 
“  20  '*  C  “"^0  I  =  1 ,  N 

IF  (A8S (X2( I ) /X( I  )-1.0 ) .G£. 1.5-7  )  GO  TO  240 
220  CONTINUE 
GC  TC  460 

240  IF  (ABS (PX1/PX2-1.  ).LE.l.E-7)  GC  TC  460 

C  FRCM^LEFTCRIGHT*  xl?  it  (  PRHV3  )  X2  (  i  )  (  PX  1  )X  (  I  )  P  X2  XlCIPl 
250  CC  260  1  =  1, N 
260  X1(I)=X(I) 

C  TFFCW  AWAY  PIGMT  PART 


F1=FX2 

C  TEMPORARILY  IN  X  STORAGE 
CALL  EVALU 

GC  TQ  (540,280,170),  NSATIS 
280  CONTINUE 

C  SWITCF^VECTCRS  TO  PROPER  POSITION 
FX1=P0 

CC  290  1=1, N 
XX  =  X2( I  ) 

X2(I)=X(I) 

X (I)=XX 


290 


GC  TG  200 


C  LEFT  'hcE 
C—  CHANGES 


TOSSED  AWAY 
FOR  NONUNI MCCAL  FN 


63 


C—  GC  7C  ThRCW  AWAY 


400 

410 

420 


420 

440 


450 


C  THE 
460 


4  70 


480 

490 

500 

C  IF 

510 
5  20 
550 

C  ARE 
540 

5i0 


C - 

c  - 


RIGHT  IN  THIS  CASEINIT  VAL  LT  FI3  PT 
350,550,410 


COMPUTE  MI 


IF  (PREV2-PX2) 

CC  420  1^1, N 

>5(ij=x2( n 

>2(I)=X(I) 

FPEV3=PX1 
F>1=PX2 
CC  440  1  =  1, N 

Xd  )=0, 28196601=)' (  XI  (  U-X2(I))+X2(I) 

CALL  eVALU 

GC  TC  (540,450,170),  NSATIS 
CONTINUE 
F>2=P0 
GC  TC  300 

INTERIOR  POINTS  NON  GIVE  EQUAL  VALUE  FOR  F, 

CC  470  1  =  1, N 
CELXOCI  )  =  X(I ) 

X (I )  =  (OELXO{I  )  +  X2(I ) )^Q.5 

CONTINUE 

CALL  EVALU 

GC  TO  (  480,490),  KSI« 

IF  (ABS(P0/PXl-l.),GT.l*E-7)  GO  TC  520 
GC  TO  (500,510),  ISW 
IF  (P0.LT.P21)  GO  TC  510 
ISW  =  2 

P-FLNCTICN  CIDN,T  GC  DOWN  TRY  NEC-  VECT. 

GC  TO  20 
RETURN 

CC  530  1  =  1, N 
X ( I )=OELXO( I ) 

GC  TO  250 

HE  NCW  IN  FEASIBILITY  PHASE 
CC  550  I=1,M 
IF  (RJ(I))  560,5o0,550 
CCNTINUE 
NSATIS=4 
RETURN 

PFCELEM  HAS  BECOME  FEASIBLE 
F  -  FUNCTION  CHANGES  IF  A  CONSTRAINT  BECCMES  FEASIBLE 


560 

MN  =  C 

CC  570 

1  =  1, M 

570 

RwKI  )  = 

PJ(I) 

RETURN 

580 

FORMAT 

(  80H 

OPT  CAN-T  FINC  A 


IGIVES  A 
ENC 


FEASIBLE  POINT, THAT 
INCTION.  ) 
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GLOSSARY  FOR  CUBIC  INTERPOLATION  OPT  SUBROUTINE 


ADELX 

AL 

BL 

COTES 

D 

DELX 

DELXO 

DLXl 

DLX2 

DOTS 

DOTT 

DOTTA 

DOTTB 

EPSO 

I 

J 

K 

KTER 

KTR 

M 

MN 


Magnitude  of  the  search  direction  vector 

Lower  step  length  0^ 

L 

Upper  step  length  0^ 

Cosine  of  the  angle  between  the  gradient  and  search 
direction  vectors 

Factor  by  which  increment  is  multiplied 
Search  direction  vector  s 

Gradient  vector  of  the  penalty  fiinction  VP(X) 

Gradient  vector  of  penalty  function  for  lower  step 
length  VP(X3) 

Gradient  vector  of  penalty  function  for  upper  step 
length  VP(X2) 

Directional  derivative  of  interpolated  step  length 
times  ADELX  P' (X) 

Directional  derivative  of  initial  starting  point 
times  ADELX  P' (XI) 

Directional  derivative  of  lower  step  length  times 
ADELX  P' (X3) 

Directional  derivative  of  upper  step  length  times 
ADELX  P' (X2) 

Stopping  criteria  tolerance  for  cosine  test 
Index  used  in  DO  loops 
Index  used  in  DO  loops 
Index  used  in  DO  loops 

Niomber  of  points  evaluated  in  bracketing  the  minimiam. 

Number  of  cubic  interpolations  performed 

Number  of  inequality  constraints 

Number  of  moves  in  search  for  solution  of  a 
subproblem 
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N 

Niomber  of  variables 

NRED 

Number  of  consecutive  step  length  reductions 

NSATIS 

Indicates  whether  constraints  are  satisfied 

NTCTR 

Number  of  the  point  on  which  the  program  is  working 

N405 

Switch  showing  that  a  feasible  point  on  the 
direction  vector  could  not  be  found  and  the 
negative  was  tried. 

PO 

Penalty  function  value  at  current  xvector  P(X) 

PXl 

Penalty  faction  value  at  lower  step  length  P{X3) 

PX2 

Penalty  function  value  at  upper  step  length  P(X2) 

Q 

Coefficient  used  in  computing  interpolated  step  length 

RJ 

Vector  of  current  values  of  the  constraints 

RJl 

Vector  of  previous  values  of  the  constraints 

SMAG 

Magnitude  of  the  gradient  of  the  penalty  function 

TAL 

Current  step  length  0 

TING 

Increment  by  which  step's  length  increased  or 
decreased 

X 

Current  X  vector 

XI 

Initial  starting  point  X  vector 

X2 

X  vector  of  upper  step  length  Xq+0^S 

X3 

X  vector  of  lower  step  length  Xq+0j^S 

Z 

Coefficient  used  in  computing  interpolated  step 
length 
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^  ^  SLEFCUTINE  OPT 

r  I  CUBIC  INTERPCLATIQN  7C  FIND  Th6 

C  ALCNG  THE  GIVEN  SEARCH  VECTCR 

INPLICIT  REAL^8(4—H  Q*"Z) 

CCHHGN/SHARE/  X(IOO),  CEL(IOO),  A ( 100, 100  )  ,N  ,  NN,NF1 

1 ,  ^  M 

?^Jil!2i^//){H^-/.-'^»C-,P0,RSIGMA,  PJ(2C0),  RHC 
CCNNGN/CRST/  DELX(IOO),  OcLXO(lOO),  RH3IN, RATIO,  EFSI, 
1 ( HETAO, 

2  PSIGl,  Gl,  Xl(lOC),  X2(100),  X3(100),  XR2(100), 
3XF1(100  ),PR1, 

i  RJKZOC),  DOTT,  PGRAC(IOC),  CIAG(IOO), 

Sr^PREVS,  ADfcLX,  NTCTR,  NUHINI,  NPHASE,  NSATIS 
CINENSICN  CLXK 100  )  ,CLX2( IOC  ) 

AES(CUHN'Y)=CABS(OUPi'<Y) 

S CRT ( DUMMY )=DSQRT( DUMMY) 

EFSC=l»E-7 
CCTS=O.Q 
N AC  5=1 
GC  TO  5 
'  DC  A  1=1, N 
A  CELX{ I) =-DeLX( I ) 

5  CONTINUE 
C  INITIALIZATION 
^^=MN+l 
NTCTR=NTCTR+1 
CCTTA=0.0 
TAL=0.0 
AL=C.O 
C  =  2.0 
KTER=0 
NFEC=0 
TINC=1,0 
SMAG=0o0 
CC  10  K=1,N 
SMAG=SMAG+CELX(K )*«2 
CCTTA  =  CCTTA-DELX(K  )=^D£LX0(K) 

X1(K)=X (K ) 

CLXK  K)  =0ELX0(K) 

1C  X'(K)=X(K) 

CCTT=-CCTTA 

FX1=P0 

C  INITIAL  BRACKETING  BEGINS 
25  TAL=TAL+TINC 
2C  CC  25  K=1,N 
25  X(K)=X1(K)+DELX(KJ«TAL 
KTER=KTER+1 
CALL  6VALU 

GC  TO  (170,50,40),  NSATIS 
C  RECLCE  STEP  SIZE 
AC  C=C.5 

NFEC=NRE0+1 
TINC=(T AL-AL)*D 
TAL=TAi.-TINC 
GC  TC  20 

C  SECCNC  FEASIBLE  POINT  FOUND 
5C  Bl=TAL 
FX2=P0 

CALL  GRA0(2) 

CCTTB=0.0 
CC  60  K=1,N 
X2(K)=X (K  ) 

CLX2(K) =D£LX0(K)  ^  ^  , 

6C  CCTTE  =  CCTTe-OELX(K )*DELXO(K  ) 

C  CHECK  STOPPING  CRITERIA  OR  MAYBE  REVERSE  SEARCH  uIRECTICN 
I f (ADELX.EC.O  .0)  GC  TO  91 
IF (SMAG.EC.0.0)  GC  TO  91 
CCTES=  ABS(DOTTB  )/(SQRT(SMAG)  =!'ACELX) 

IF  (COTES.LTaEPSQ  )  GO  TO  91 
IF (KTER.GT.lOO)  GC,TC^280 
I F (NREC  .GT.20)  GO  TO  280 
CC  75  K=1,N 
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7*=  ^'’ 

GC  TO  260 

76^IF(4eS(ABS(PX2/PXl)-l. ).L£.l.£-7)  GO  TO  91 
C  If  CEPIVATIVE  IS  POSITIVE  THEN  ThE  MINIMUM  IS  5FACK5TEC 
C  ANC  CLEIC  INTERPOLATION  CAN  BE  PERFORMED,  IF  NGT  MAKE  THIS 
C  PCINT  TFE  LOWER  STEP  LENGTH  AND  STEF  OUT  FARTHER 
IF(DOTTE«GT.O.O)  GO  TO  100 
N  R  E  C  =  0 
TINC=TINC*D 
CO  60  K=1,N 
DL>1(K)=0LX2(K) 

6C  X2(K)=X2(K) 

F>1=PX2 
AL  =  EL 

CCTTA=COTTB 
GO  TO  25 

C  OF  THE  TWO  FEASIBLE  POINTS  RETURN  WITH  THE  SMALLEST 

51  IF(FX2.LE.FX1)  GC  TO  220 
FC=FX1 

CO  S2  K=1,N 
X (K )  =  X3  (K  ) 

CELX0(K)=DLX1(K) 

52  CONTINUE 
GC  TO  220 

C  OLEIC  interpolation  FOLLOWS 
ICO  KTF=0 
105  CONTINUE 

Z  =  3.=i‘((  FX1-PX2)/(EL-AL)  ) +D0TTA  +  CCTT8 
C=SCRT(Z**2-DOTTA*CCTT8 ) 

TAL  =  BL-  (BL-AD^IDCTTB  +  Q-Z)/  (D0TTB-CGTTA  +  2.=»C  ) 

CO  130  l<  =  l,N 

12C  X  (K)  =  X1  (K)+OELX(K)=i'TAL 
CALL  EVALU 

GC  TO  (170,135,40),  NSATIS 
135  CALL  GRA0(2) 

CCTS=O.Q 
CC  140  K=1,N 

CCTS=CCTS-CELX(K)*0£LX0(K) 
lAC  CONTINUE 

C  CHECK  STOPPING  CRITERION  CR  IF  SEARCH  DIRECTION  SHOULD  EE 
C  REVERSEC 

I F(ADELX.EC.O.O)  GO  TO  471 
IF(SMAG.Ee.O.O)  GC  TO  471 
CCTES=AES( COTS)/ (SCR7(SMAG) «AD£LX ) 

IF  (COTES. LT. EPSQ )  GO  TO  471 

IF(ABS(ABS  (PX2/P0  ) -1.  )  .L  E.  l<.E-7  )  GO  TO  471 

KTR=KTR+1 

IF  (KTR.GT.IO)  GO  TO  471 
CC  143  K  =  1,N 

IF(AeS(A8S(X2(K)/X(K)  ) -1 . ) .GT . 1 , E-7 )  GO  TC  163 
143  CCNTINUE 
GC  TC  220 

163  IF(COTS.GT.O.O)  GO  TO  150 
C  THROW  AWAY  LEFT  SIDE 
Al=TAL 
FX1=PC 
CCTTA=CCTS 
CC  145  K=1,N 
CLX1(K)=DELX0(K  ) 

145  X3(K)=X(K) 

GC  TO  105 

C  THFCW  AWAY  RIGHT  SIDE 
15C  eL=TAL 
FX2*P0 
CCTT8=CCTS 
CC  160  K=1,N^^  ^ 

CLX2(K)=D£LX0(K) 

16C  X2(K}=X(K) 

GC  TO  105 
170  CC  180  K=1,M 

IF(RJ(K))  190,190,160 
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16C  CCNTINtE 
NS^TIS=4 
RETURN 
ISC  NN=0 

CC  200  K  =  1,M 
2C0  RJ1(K)=PJ(K) 

RETURN 

RETURN  l^ITH  SMAUUEST  VAUUc 

471  I  F { FO.UE.PXl )  GO  TO  474 
IF(FX1.UE.PX2)  GO  TO  475 

472  FC=FX2 

CC  473  K=l,N 
X(K)=X2(K) 

CEUXO(K  )  =  DUX2(K) 

473  CCNTINUE 
GC  TO  220 

474  IF(F0.UE«PX2)  GO  TC  220 
GC  TC  472 

475  PC=FX1 

CC  476  K=1,N 
X  (K  )  =  X3 (K) 

CEUXO<K  )=DUX1(K) 

476  CCNTINUE 
22C  CCNTINUE 
225  CCNTINUE 

RETURN 

REVERSE  SEARCH  OIRECTICNS 
26C  GC  TO  ( 290,200) ,N405 
29C  N4C5=2 

NTCTR=NTCTR-1 
NN=NN-l 
GC  TO  3 

3CC  NFITE(6,580) 

CALU  TIHEC 
C6UU  QUTPU'T(l) 

CAUU  REJECT 

560  FCFMAT(  80H  OPT  CAN-T  FIND  A  FEASIfiUE  PCINT,THAT 
IGIVES  A  UOWER  VAUUE  OF  THE  F-FUNCTIQN.  ) 

STCF  22042 
ENC 
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GLOSSARY  FOR  QUADRATIC  INTERPOLATION  OPT  SUBROUTINE 


AL 

Interpolated  step  length 

D 

Factor  by  which  increment  is  multiplied 

DELX 

Vector  indicating  the  direction  of  move  in  one 
dimensional  optimization  S 

DELXO 

Gradient  vector  of  the  penalty  function  VP 

DOTT 

Directional  derivative  of  initial  starting  point 
times  search  direction  magnitude  P'(Xl) 

I 

Index  used  in  DO  loops 

J 

Index  used  in  DO  loops 

K 

Index  used  in  DO  loops 

KTER 

Nxxmber  of  points  evaluated  in  bracketing  the 
minimum 

KTR 

N\amber  of  quadratic  interpolations  performed 

M 

Nximber  of  inequality  constraints 

MN 

Number  of  moves  in  search  for  solution  of  a 
subproblem 

N 

Number  of  variables 

NFIRS 

Indicates  if  more  than  two  feasible  step  lengths 
found 

NRED 

Number  of  consecutive  step  length  reductions 

NS ATI S 

Indicates  whether  constraints  are  satisfied 

NTCTR 

Number  of  the  point  on  which  the  program  is  working 

N405 

Switch  showing  that  a  feasible  point  on  the  direction 
vector  could  not  be  found  and  the  negative  was  tried 

PO 

Current  penalty  foonction  value  P{X) 

PXl 

Penalty  fiinction  value  of  lower  step  length  P(X3) 
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PX2 

RJ 

TAL 

TALI 

TAL  2 

TINC 

X 

XX 

XI 
X2 
X3 


Penalty  fvmction  value  of  upper  step  length 
Vector  of  previous  values  of  the  constraints 
Current  step  length  0 
Lower  step  length  0^. 

Li 

Upper  step  length  0^ 

Increment  by  which  step  length  increased  or 
decreased 

Current  X  vector  Xq+  0s 

Temporary  storage  used  for  switching  values 
Initial  starting  point  vector  Xq 
X  vector  at  upper  step  length  ^o^®u^ 

X  vector  of  lower  step  length  Xq+0j^s 
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4 

c 


SIEROUTINE  OPT 

C  This  CPT  SUEROUTINS  FINDS  THE  MINIMUM  ALONG  A  GIVEN  SEAPCH 
C  VeCTCF  USING  CUAORATIC  INTERPOLATION 
IMPLICIT  SEAL’^'Si  A-H,a-Z  ) 

CCMMON/SHARE/  X(IOO),  DELdOO),  A  ( 100 1 100  )  ,  N  ,  V  ,  PN,NF1 
ll^^'l 

LCPNCN  /VALUE/  F , G  ,  PO , RS I GM A ,  RJ(200)»  RFC 
CCMMGN/CRST/  DELX(100)t  0ELXC(100)»  RHQIN, PATIO,  EPSI, 
ITFETAO, 

2  PSIGl,  Gl,  XKIOO),  X2(100),  >3(100),  X02(iOC), 
EXFKIQO  ),PR1, 

4  FF2,F1,  FI,  RJ1(200),  DOTT,  PGRAC(IOO),  CIAGdOO), 

5  FP£V3,A0ELX,  NTCTR,  NUKINI,  NPHASE,  NSATIS 
AES  (CUN'MY)=DABS(DLN'NY) 

N405=l 

CCTT=C. 

C  C  2  K  1  N 

2  CCTT  =  CCTf+DELX(K  )=^CELXO(K) 

GO  TO  5 
2  CC  4  1=1, N 

C£LXd)=-DELX(I  ) 

CONTINUE 
C  INITIALIZATION 
NN=NN+1 
NTCTR=NTCTR+l 
TAL=0.0 
TAL1=0.0 
C  =  2. 

KTER=0 
NREC=0 
TINC=1.0 
CC  7  K*1,N 
X1(K)=X(K) 

X3(K)=X (K) 

7  FX1=F0 

NFIRS=0 

C  START  ERACKETING  PROCEDURE 
TAL=TAL4TINC 
10  CC  13  K=1,N 

12  X (K  )  =  X1 (K)+OELX(K  )*TAL 

KT£P=KTER+1 
CALL  EVALU 

GC  TO  ( 170,20,20) ,  NSATIS 
C  RECUCE  STEP  SIZE 
20  C=C.5 

IF  (NREC.GT.20)  GO  TO  280 
NREC=NRED+1 
TINC=(TAL-TAL1)#0 
T  AL  =  TAL-TINC 
GC  TO  1C 

C  HAVE  2  CR  MCRE  FEASIBLE  POINTS  BEEN  FCUND 

20  IFiNFIRS.EC.l)  GO  TC  40 
NFIRS=1 

21  CC  32  K=1,N 

22  X2(K)=X(Kj 
FX2=P0 

IF(FX2.G£.PX1)  GO  TO  20 
C  INCREASE  STEP  SIZE 
25  TINC=TINC*D 

TAL=TAL+TINC 
NF£C=0 
GC  TO  10 
40  CONTINUE 

C  RECRCER  POINTS  „ 

IF(TAL. LT.TAL2)  GO  TO  54 

IF(F0«LT.PX2)  60  TO  50 

X>=TAL 

TAL=TAL2 

TAL2=XX 

X>  =  FO 

FC=FX2 
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F>2=XX 

CC  45  K=1,N 

XX=X(K) 

X(K)=X2{K) 

45 

X2{K)=XX 

GC  TO  54 

46 

TAL1=TAL 

P>1=P0 

CC  47  K=1,N 
X3 (K)=X (K) 
TAL=TAL2 

GC  TO  35 

47 

50 

TAL1=TAL2 

F>1=PX2 

CC  51  K=1,N 
X3(K)=X2(K ) 

51 

>2(K)=X(K) 

TAL2=TAL 

FX2=P0 

52 

B  -a 

c 

c 

54 

c  c 

56 


C  CHECK  STOPPING  CRITERION  OR  IF  SEARCH  DIRECTION  SHGULC  EE 
C  REVERSEC 

IF(AeS(  ABS{  PO/PXl  »-l.  )  .LT.l.e-7  )  C-C  TO  220 
CC  52  K=1.N 

IF(ABS{ ABS(X1(K)/X{K)  )-l. ).GT.l.E-7)  GO  TC  53 
CCNTINLE 
GC  TO  260 
CCNTINLE 
GC  TO  35 

CHECK  STOPPING  CRITERION  OR  IF  SEARCH  CIRECTION  SHQULC  BE 
REVERSED 

IF{AeS{ABS(P0/PXl)-l.).LT.l.E-7)  GO  TO  220 
CC  55  K=1,N 

IF(ABS{ABS{Xl(KJ/X{Kn-l.).GT.l.E-7)  GO  TC  56 
CONTINUE 
GC  TO  260 
CCNTINLE 

IF  (KTER.GE.lOO)  GO  TO  280 
IF(F0.GT.PX2)  GO  TC  46 
IF(FO.GE.PXl)  GO  TO  31 
CC  57  K=1,N 

IF(ABS(ABS(X2{K)/X{K)  )-l.).G6.1.E-7)  GO  TO  60 
57  CONTINUE 

GC  TC  220 

C  CLADRATIC  INTERPOLATION 

60  KTR=0 

61  KTP— 

AL=0.5«  ( (  TAL**2-TAL2=f‘*2)*PXl+  (TA  L  2*’*‘2-T  AL  1=)'=>2  )  »P0  + 

1  (TALl^^E-TAL^’Ka 

2  )  ♦FX2)  /  { (TAL-TAL2  )=«‘PX1+  (T AL 2-TAL 1 )  =S'PO+ (  T AL  1-T AL  ) P X2  ) 
IFdAL.GT.AL)  GO  TC  75 

C  THPCU  AWAY  LEFT  SIDE 
TAL1=TAL 
FX1=P0 
CC  71  K=1,N 
71  X3(K)=X(K) 

GC  TO  78 

C  T|.pcw  AWAY  RIGHT  SICt 
75  TAL2=TAL 

FX2=P0 
CC  77  K=1,N 

77  X2(K)=X(K) 

78  TAL=AL 

CC  79  K=lrN 

79  X (K )  =  X1 (K)+OELX(K  )*TAL 
CALL  EVALU 

GC  TC  (  170,80 »85)  ,  NSATIS 

C°CH6CK*’ST0PPINd^CRITERI0N  OR  IF  SEARCH  DIRECTION  SHGULC  8  = 
r  rpC  =  r<!FC 

IF{AES{ABS{X2{K)/X(K) )-l. ).GE.l.E-7)  GO  TC  53 
82  CCNTINUE 

GC  TO  220 
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IF  5|S(4BS(PO/PX1)-1o).LT,1.E-7)  goto  220 
IF(KTR.GT.20)  GO  TO  220 
GC  TO  61 
65  NFEC=0 
KTtB=0 
GC  TO  20 

170  CC  180  K=l,y 

IF(BJ(KJ)  190,190,180 
160  CONTINUE 
NS^TIS=4 
RETURN 
190  MN=0 

CC  200  K=1,M 
200  RJ1(K)=PJ(K) 

RETURN 

220  CONTINUE 

C  RETURN  NITH  SMALLEST  VALUE 
IFIFXl.GE.PO)  GO  TO  2A0 
PC=FX1 

CC  230  K=1,N 
220  X(K)=X3(K» 

2A0  CONTINUE 
R  ET UR N 

C  REVERSE  SEARCH  DIRECTIONS 
26C  GC  TO  ( 290,300) ,N405 
290  N405=2 

NTCTR=NTCTR-1 
NN=MN-1 
GC  TG  3 

300  WPITE{6,580) 

CALL  TINEC 
CALL  OUTPUT(i) 

CALL  REJECT 

560  FGFMAT(  80H  OPT  CAN-T  FIND  A  FEASIBLE  PCINT,ThAT 
IGIVES  A  LOWER  VALUE  OF  THE  P-FUNCTIGN.  ) 

STOP  22042 
ENC 
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